home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / zlatrs.f < prev    next >
Text File  |  1996-07-19  |  28KB  |  881 lines
  1.       SUBROUTINE ZLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
  2.      $                   CNORM, INFO )
  3. *
  4. *  -- LAPACK auxiliary routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     June 30, 1992
  8. *
  9. *     .. Scalar Arguments ..
  10.       CHARACTER          DIAG, NORMIN, TRANS, UPLO
  11.       INTEGER            INFO, LDA, N
  12.       DOUBLE PRECISION   SCALE
  13. *     ..
  14. *     .. Array Arguments ..
  15.       DOUBLE PRECISION   CNORM( * )
  16.       COMPLEX*16         A( LDA, * ), X( * )
  17. *     ..
  18. *
  19. *  Purpose
  20. *  =======
  21. *
  22. *  ZLATRS solves one of the triangular systems
  23. *
  24. *     A * x = s*b,  A**T * x = s*b,  or  A**H * x = s*b,
  25. *
  26. *  with scaling to prevent overflow.  Here A is an upper or lower
  27. *  triangular matrix, A**T denotes the transpose of A, A**H denotes the
  28. *  conjugate transpose of A, x and b are n-element vectors, and s is a
  29. *  scaling factor, usually less than or equal to 1, chosen so that the
  30. *  components of x will be less than the overflow threshold.  If the
  31. *  unscaled problem will not cause overflow, the Level 2 BLAS routine
  32. *  ZTRSV is called. If the matrix A is singular (A(j,j) = 0 for some j),
  33. *  then s is set to 0 and a non-trivial solution to A*x = 0 is returned.
  34. *
  35. *  Arguments
  36. *  =========
  37. *
  38. *  UPLO    (input) CHARACTER*1
  39. *          Specifies whether the matrix A is upper or lower triangular.
  40. *          = 'U':  Upper triangular
  41. *          = 'L':  Lower triangular
  42. *
  43. *  TRANS   (input) CHARACTER*1
  44. *          Specifies the operation applied to A.
  45. *          = 'N':  Solve A * x = s*b     (No transpose)
  46. *          = 'T':  Solve A**T * x = s*b  (Transpose)
  47. *          = 'C':  Solve A**H * x = s*b  (Conjugate transpose)
  48. *
  49. *  DIAG    (input) CHARACTER*1
  50. *          Specifies whether or not the matrix A is unit triangular.
  51. *          = 'N':  Non-unit triangular
  52. *          = 'U':  Unit triangular
  53. *
  54. *  NORMIN  (input) CHARACTER*1
  55. *          Specifies whether CNORM has been set or not.
  56. *          = 'Y':  CNORM contains the column norms on entry
  57. *          = 'N':  CNORM is not set on entry.  On exit, the norms will
  58. *                  be computed and stored in CNORM.
  59. *
  60. *  N       (input) INTEGER
  61. *          The order of the matrix A.  N >= 0.
  62. *
  63. *  A       (input) COMPLEX*16 array, dimension (LDA,N)
  64. *          The triangular matrix A.  If UPLO = 'U', the leading n by n
  65. *          upper triangular part of the array A contains the upper
  66. *          triangular matrix, and the strictly lower triangular part of
  67. *          A is not referenced.  If UPLO = 'L', the leading n by n lower
  68. *          triangular part of the array A contains the lower triangular
  69. *          matrix, and the strictly upper triangular part of A is not
  70. *          referenced.  If DIAG = 'U', the diagonal elements of A are
  71. *          also not referenced and are assumed to be 1.
  72. *
  73. *  LDA     (input) INTEGER
  74. *          The leading dimension of the array A.  LDA >= max (1,N).
  75. *
  76. *  X       (input/output) COMPLEX*16 array, dimension (N)
  77. *          On entry, the right hand side b of the triangular system.
  78. *          On exit, X is overwritten by the solution vector x.
  79. *
  80. *  SCALE   (output) DOUBLE PRECISION
  81. *          The scaling factor s for the triangular system
  82. *             A * x = s*b,  A**T * x = s*b,  or  A**H * x = s*b.
  83. *          If SCALE = 0, the matrix A is singular or badly scaled, and
  84. *          the vector x is an exact or approximate solution to A*x = 0.
  85. *
  86. *  CNORM   (input or output) DOUBLE PRECISION array, dimension (N)
  87. *
  88. *          If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
  89. *          contains the norm of the off-diagonal part of the j-th column
  90. *          of A.  If TRANS = 'N', CNORM(j) must be greater than or equal
  91. *          to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
  92. *          must be greater than or equal to the 1-norm.
  93. *
  94. *          If NORMIN = 'N', CNORM is an output argument and CNORM(j)
  95. *          returns the 1-norm of the offdiagonal part of the j-th column
  96. *          of A.
  97. *
  98. *  INFO    (output) INTEGER
  99. *          = 0:  successful exit
  100. *          < 0:  if INFO = -k, the k-th argument had an illegal value
  101. *
  102. *  Further Details
  103. *  ======= =======
  104. *
  105. *  A rough bound on x is computed; if that is less than overflow, ZTRSV
  106. *  is called, otherwise, specific code is used which checks for possible
  107. *  overflow or divide-by-zero at every operation.
  108. *
  109. *  A columnwise scheme is used for solving A*x = b.  The basic algorithm
  110. *  if A is lower triangular is
  111. *
  112. *       x[1:n] := b[1:n]
  113. *       for j = 1, ..., n
  114. *            x(j) := x(j) / A(j,j)
  115. *            x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
  116. *       end
  117. *
  118. *  Define bounds on the components of x after j iterations of the loop:
  119. *     M(j) = bound on x[1:j]
  120. *     G(j) = bound on x[j+1:n]
  121. *  Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
  122. *
  123. *  Then for iteration j+1 we have
  124. *     M(j+1) <= G(j) / | A(j+1,j+1) |
  125. *     G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
  126. *            <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
  127. *
  128. *  where CNORM(j+1) is greater than or equal to the infinity-norm of
  129. *  column j+1 of A, not counting the diagonal.  Hence
  130. *
  131. *     G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
  132. *                  1<=i<=j
  133. *  and
  134. *
  135. *     |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
  136. *                                   1<=i< j
  137. *
  138. *  Since |x(j)| <= M(j), we use the Level 2 BLAS routine ZTRSV if the
  139. *  reciprocal of the largest M(j), j=1,..,n, is larger than
  140. *  max(underflow, 1/overflow).
  141. *
  142. *  The bound on x(j) is also used to determine when a step in the
  143. *  columnwise method can be performed without fear of overflow.  If
  144. *  the computed bound is greater than a large constant, x is scaled to
  145. *  prevent overflow, but if the bound overflows, x is set to 0, x(j) to
  146. *  1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
  147. *
  148. *  Similarly, a row-wise scheme is used to solve A**T *x = b  or
  149. *  A**H *x = b.  The basic algorithm for A upper triangular is
  150. *
  151. *       for j = 1, ..., n
  152. *            x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
  153. *       end
  154. *
  155. *  We simultaneously compute two bounds
  156. *       G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
  157. *       M(j) = bound on x(i), 1<=i<=j
  158. *
  159. *  The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
  160. *  add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
  161. *  Then the bound on x(j) is
  162. *
  163. *       M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
  164. *
  165. *            <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
  166. *                      1<=i<=j
  167. *
  168. *  and we can safely call ZTRSV if 1/M(n) and 1/G(n) are both greater
  169. *  than max(underflow, 1/overflow).
  170. *
  171. *  =====================================================================
  172. *
  173. *     .. Parameters ..
  174.       DOUBLE PRECISION   ZERO, HALF, ONE, TWO
  175.       PARAMETER          ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0,
  176.      $                   TWO = 2.0D+0 )
  177. *     ..
  178. *     .. Local Scalars ..
  179.       LOGICAL            NOTRAN, NOUNIT, UPPER
  180.       INTEGER            I, IMAX, J, JFIRST, JINC, JLAST
  181.       DOUBLE PRECISION   BIGNUM, GROW, REC, SMLNUM, TJJ, TMAX, TSCAL,
  182.      $                   XBND, XJ, XMAX
  183.       COMPLEX*16         CSUMJ, TJJS, USCAL, ZDUM
  184. *     ..
  185. *     .. External Functions ..
  186.       LOGICAL            LSAME
  187.       INTEGER            IDAMAX, IZAMAX
  188.       DOUBLE PRECISION   DLAMCH, DZASUM
  189.       COMPLEX*16         ZDOTC, ZDOTU, ZLADIV
  190.       EXTERNAL           LSAME, IDAMAX, IZAMAX, DLAMCH, DZASUM, ZDOTC,
  191.      $                   ZDOTU, ZLADIV
  192. *     ..
  193. *     .. External Subroutines ..
  194.       EXTERNAL           DLABAD, DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTRSV
  195. *     ..
  196. *     .. Intrinsic Functions ..
  197.       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
  198. *     ..
  199. *     .. Statement Functions ..
  200.       DOUBLE PRECISION   CABS1, CABS2
  201. *     ..
  202. *     .. Statement Function definitions ..
  203.       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
  204.       CABS2( ZDUM ) = ABS( DBLE( ZDUM ) / 2.D0 ) +
  205.      $                ABS( DIMAG( ZDUM ) / 2.D0 )
  206. *     ..
  207. *     .. Executable Statements ..
  208. *
  209.       INFO = 0
  210.       UPPER = LSAME( UPLO, 'U' )
  211.       NOTRAN = LSAME( TRANS, 'N' )
  212.       NOUNIT = LSAME( DIAG, 'N' )
  213. *
  214. *     Test the input parameters.
  215. *
  216.       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  217.          INFO = -1
  218.       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
  219.      $         LSAME( TRANS, 'C' ) ) THEN
  220.          INFO = -2
  221.       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
  222.          INFO = -3
  223.       ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT.
  224.      $         LSAME( NORMIN, 'N' ) ) THEN
  225.          INFO = -4
  226.       ELSE IF( N.LT.0 ) THEN
  227.          INFO = -5
  228.       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  229.          INFO = -7
  230.       END IF
  231.       IF( INFO.NE.0 ) THEN
  232.          CALL XERBLA( 'ZLATRS', -INFO )
  233.          RETURN
  234.       END IF
  235. *
  236. *     Quick return if possible
  237. *
  238.       IF( N.EQ.0 )
  239.      $   RETURN
  240. *
  241. *     Determine machine dependent parameters to control overflow.
  242. *
  243.       SMLNUM = DLAMCH( 'Safe minimum' )
  244.       BIGNUM = ONE / SMLNUM
  245.       CALL DLABAD( SMLNUM, BIGNUM )
  246.       SMLNUM = SMLNUM / DLAMCH( 'Precision' )
  247.       BIGNUM = ONE / SMLNUM
  248.       SCALE = ONE
  249. *
  250.       IF( LSAME( NORMIN, 'N' ) ) THEN
  251. *
  252. *        Compute the 1-norm of each column, not including the diagonal.
  253. *
  254.          IF( UPPER ) THEN
  255. *
  256. *           A is upper triangular.
  257. *
  258.             DO 10 J = 1, N
  259.                CNORM( J ) = DZASUM( J-1, A( 1, J ), 1 )
  260.    10       CONTINUE
  261.          ELSE
  262. *
  263. *           A is lower triangular.
  264. *
  265.             DO 20 J = 1, N - 1
  266.                CNORM( J ) = DZASUM( N-J, A( J+1, J ), 1 )
  267.    20       CONTINUE
  268.             CNORM( N ) = ZERO
  269.          END IF
  270.       END IF
  271. *
  272. *     Scale the column norms by TSCAL if the maximum element in CNORM is
  273. *     greater than BIGNUM/2.
  274. *
  275.       IMAX = IDAMAX( N, CNORM, 1 )
  276.       TMAX = CNORM( IMAX )
  277.       IF( TMAX.LE.BIGNUM*HALF ) THEN
  278.          TSCAL = ONE
  279.       ELSE
  280.          TSCAL = HALF / ( SMLNUM*TMAX )
  281.          CALL DSCAL( N, TSCAL, CNORM, 1 )
  282.       END IF
  283. *
  284. *     Compute a bound on the computed solution vector to see if the
  285. *     Level 2 BLAS routine ZTRSV can be used.
  286. *
  287.       XMAX = ZERO
  288.       DO 30 J = 1, N
  289.          XMAX = MAX( XMAX, CABS2( X( J ) ) )
  290.    30 CONTINUE
  291.       XBND = XMAX
  292. *
  293.       IF( NOTRAN ) THEN
  294. *
  295. *        Compute the growth in A * x = b.
  296. *
  297.          IF( UPPER ) THEN
  298.             JFIRST = N
  299.             JLAST = 1
  300.             JINC = -1
  301.          ELSE
  302.             JFIRST = 1
  303.             JLAST = N
  304.             JINC = 1
  305.          END IF
  306. *
  307.          IF( TSCAL.NE.ONE ) THEN
  308.             GROW = ZERO
  309.             GO TO 60
  310.          END IF
  311. *
  312.          IF( NOUNIT ) THEN
  313. *
  314. *           A is non-unit triangular.
  315. *
  316. *           Compute GROW = 1/G(j) and XBND = 1/M(j).
  317. *           Initially, G(0) = max{x(i), i=1,...,n}.
  318. *
  319.             GROW = HALF / MAX( XBND, SMLNUM )
  320.             XBND = GROW
  321.             DO 40 J = JFIRST, JLAST, JINC
  322. *
  323. *              Exit the loop if the growth factor is too small.
  324. *
  325.                IF( GROW.LE.SMLNUM )
  326.      $            GO TO 60
  327. *
  328.                TJJS = A( J, J )
  329.                TJJ = CABS1( TJJS )
  330. *
  331.                IF( TJJ.GE.SMLNUM ) THEN
  332. *
  333. *                 M(j) = G(j-1) / abs(A(j,j))
  334. *
  335.                   XBND = MIN( XBND, MIN( ONE, TJJ )*GROW )
  336.                ELSE
  337. *
  338. *                 M(j) could overflow, set XBND to 0.
  339. *
  340.                   XBND = ZERO
  341.                END IF
  342. *
  343.                IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN
  344. *
  345. *                 G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
  346. *
  347.                   GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) )
  348.                ELSE
  349. *
  350. *                 G(j) could overflow, set GROW to 0.
  351. *
  352.                   GROW = ZERO
  353.                END IF
  354.    40       CONTINUE
  355.             GROW = XBND
  356.          ELSE
  357. *
  358. *           A is unit triangular.
  359. *
  360. *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
  361. *
  362.             GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
  363.             DO 50 J = JFIRST, JLAST, JINC
  364. *
  365. *              Exit the loop if the growth factor is too small.
  366. *
  367.                IF( GROW.LE.SMLNUM )
  368.      $            GO TO 60
  369. *
  370. *              G(j) = G(j-1)*( 1 + CNORM(j) )
  371. *
  372.                GROW = GROW*( ONE / ( ONE+CNORM( J ) ) )
  373.    50       CONTINUE
  374.          END IF
  375.    60    CONTINUE
  376. *
  377.       ELSE
  378. *
  379. *        Compute the growth in A**T * x = b  or  A**H * x = b.
  380. *
  381.          IF( UPPER ) THEN
  382.             JFIRST = 1
  383.             JLAST = N
  384.             JINC = 1
  385.          ELSE
  386.             JFIRST = N
  387.             JLAST = 1
  388.             JINC = -1
  389.          END IF
  390. *
  391.          IF( TSCAL.NE.ONE ) THEN
  392.             GROW = ZERO
  393.             GO TO 90
  394.          END IF
  395. *
  396.          IF( NOUNIT ) THEN
  397. *
  398. *           A is non-unit triangular.
  399. *
  400. *           Compute GROW = 1/G(j) and XBND = 1/M(j).
  401. *           Initially, M(0) = max{x(i), i=1,...,n}.
  402. *
  403.             GROW = HALF / MAX( XBND, SMLNUM )
  404.             XBND = GROW
  405.             DO 70 J = JFIRST, JLAST, JINC
  406. *
  407. *              Exit the loop if the growth factor is too small.
  408. *
  409.                IF( GROW.LE.SMLNUM )
  410.      $            GO TO 90
  411. *
  412. *              G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
  413. *
  414.                XJ = ONE + CNORM( J )
  415.                GROW = MIN( GROW, XBND / XJ )
  416. *
  417.                TJJS = A( J, J )
  418.                TJJ = CABS1( TJJS )
  419. *
  420.                IF( TJJ.GE.SMLNUM ) THEN
  421. *
  422. *                 M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
  423. *
  424.                   IF( XJ.GT.TJJ )
  425.      $               XBND = XBND*( TJJ / XJ )
  426.                ELSE
  427. *
  428. *                 M(j) could overflow, set XBND to 0.
  429. *
  430.                   XBND = ZERO
  431.                END IF
  432.    70       CONTINUE
  433.             GROW = MIN( GROW, XBND )
  434.          ELSE
  435. *
  436. *           A is unit triangular.
  437. *
  438. *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
  439. *
  440.             GROW = MIN( ONE, HALF / MAX( XBND, SMLNUM ) )
  441.             DO 80 J = JFIRST, JLAST, JINC
  442. *
  443. *              Exit the loop if the growth factor is too small.
  444. *
  445.                IF( GROW.LE.SMLNUM )
  446.      $            GO TO 90
  447. *
  448. *              G(j) = ( 1 + CNORM(j) )*G(j-1)
  449. *
  450.                XJ = ONE + CNORM( J )
  451.                GROW = GROW / XJ
  452.    80       CONTINUE
  453.          END IF
  454.    90    CONTINUE
  455.       END IF
  456. *
  457.       IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN
  458. *
  459. *        Use the Level 2 BLAS solve if the reciprocal of the bound on
  460. *        elements of X is not too small.
  461. *
  462.          CALL ZTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
  463.       ELSE
  464. *
  465. *        Use a Level 1 BLAS solve, scaling intermediate results.
  466. *
  467.          IF( XMAX.GT.BIGNUM*HALF ) THEN
  468. *
  469. *           Scale X so that its components are less than or equal to
  470. *           BIGNUM in absolute value.
  471. *
  472.             SCALE = ( BIGNUM*HALF ) / XMAX
  473.             CALL ZDSCAL( N, SCALE, X, 1 )
  474.             XMAX = BIGNUM
  475.          ELSE
  476.             XMAX = XMAX*TWO
  477.          END IF
  478. *
  479.          IF( NOTRAN ) THEN
  480. *
  481. *           Solve A * x = b
  482. *
  483.             DO 120 J = JFIRST, JLAST, JINC
  484. *
  485. *              Compute x(j) = b(j) / A(j,j), scaling x if necessary.
  486. *
  487.                XJ = CABS1( X( J ) )
  488.                IF( NOUNIT ) THEN
  489.                   TJJS = A( J, J )*TSCAL
  490.                ELSE
  491.                   TJJS = TSCAL
  492.                   IF( TSCAL.EQ.ONE )
  493.      $               GO TO 110
  494.                END IF
  495.                TJJ = CABS1( TJJS )
  496.                IF( TJJ.GT.SMLNUM ) THEN
  497. *
  498. *                    abs(A(j,j)) > SMLNUM:
  499. *
  500.                   IF( TJJ.LT.ONE ) THEN
  501.                      IF( XJ.GT.TJJ*BIGNUM ) THEN
  502. *
  503. *                          Scale x by 1/b(j).
  504. *
  505.                         REC = ONE / XJ
  506.                         CALL ZDSCAL( N, REC, X, 1 )
  507.                         SCALE = SCALE*REC
  508.                         XMAX = XMAX*REC
  509.                      END IF
  510.                   END IF
  511.                   X( J ) = ZLADIV( X( J ), TJJS )
  512.                   XJ = CABS1( X( J ) )
  513.                ELSE IF( TJJ.GT.ZERO ) THEN
  514. *
  515. *                    0 < abs(A(j,j)) <= SMLNUM:
  516. *
  517.                   IF( XJ.GT.TJJ*BIGNUM ) THEN
  518. *
  519. *                       Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM
  520. *                       to avoid overflow when dividing by A(j,j).
  521. *
  522.                      REC = ( TJJ*BIGNUM ) / XJ
  523.                      IF( CNORM( J ).GT.ONE ) THEN
  524. *
  525. *                          Scale by 1/CNORM(j) to avoid overflow when
  526. *                          multiplying x(j) times column j.
  527. *
  528.                         REC = REC / CNORM( J )
  529.                      END IF
  530.                      CALL ZDSCAL( N, REC, X, 1 )
  531.                      SCALE = SCALE*REC
  532.                      XMAX = XMAX*REC
  533.                   END IF
  534.                   X( J ) = ZLADIV( X( J ), TJJS )
  535.                   XJ = CABS1( X( J ) )
  536.                ELSE
  537. *
  538. *                    A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
  539. *                    scale = 0, and compute a solution to A*x = 0.
  540. *
  541.                   DO 100 I = 1, N
  542.                      X( I ) = ZERO
  543.   100             CONTINUE
  544.                   X( J ) = ONE
  545.                   XJ = ONE
  546.                   SCALE = ZERO
  547.                   XMAX = ZERO
  548.                END IF
  549.   110          CONTINUE
  550. *
  551. *              Scale x if necessary to avoid overflow when adding a
  552. *              multiple of column j of A.
  553. *
  554.                IF( XJ.GT.ONE ) THEN
  555.                   REC = ONE / XJ
  556.                   IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN
  557. *
  558. *                    Scale x by 1/(2*abs(x(j))).
  559. *
  560.                      REC = REC*HALF
  561.                      CALL ZDSCAL( N, REC, X, 1 )
  562.                      SCALE = SCALE*REC
  563.                   END IF
  564.                ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN
  565. *
  566. *                 Scale x by 1/2.
  567. *
  568.                   CALL ZDSCAL( N, HALF, X, 1 )
  569.                   SCALE = SCALE*HALF
  570.                END IF
  571. *
  572.                IF( UPPER ) THEN
  573.                   IF( J.GT.1 ) THEN
  574. *
  575. *                    Compute the update
  576. *                       x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
  577. *
  578.                      CALL ZAXPY( J-1, -X( J )*TSCAL, A( 1, J ), 1, X,
  579.      $                           1 )
  580.                      I = IZAMAX( J-1, X, 1 )
  581.                      XMAX = CABS1( X( I ) )
  582.                   END IF
  583.                ELSE
  584.                   IF( J.LT.N ) THEN
  585. *
  586. *                    Compute the update
  587. *                       x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
  588. *
  589.                      CALL ZAXPY( N-J, -X( J )*TSCAL, A( J+1, J ), 1,
  590.      $                           X( J+1 ), 1 )
  591.                      I = J + IZAMAX( N-J, X( J+1 ), 1 )
  592.                      XMAX = CABS1( X( I ) )
  593.                   END IF
  594.                END IF
  595.   120       CONTINUE
  596. *
  597.          ELSE IF( LSAME( TRANS, 'T' ) ) THEN
  598. *
  599. *           Solve A**T * x = b
  600. *
  601.             DO 170 J = JFIRST, JLAST, JINC
  602. *
  603. *              Compute x(j) = b(j) - sum A(k,j)*x(k).
  604. *                                    k<>j
  605. *
  606.                XJ = CABS1( X( J ) )
  607.                USCAL = TSCAL
  608.                REC = ONE / MAX( XMAX, ONE )
  609.                IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
  610. *
  611. *                 If x(j) could overflow, scale x by 1/(2*XMAX).
  612. *
  613.                   REC = REC*HALF
  614.                   IF( NOUNIT ) THEN
  615.                      TJJS = A( J, J )*TSCAL
  616.                   ELSE
  617.                      TJJS = TSCAL
  618.                   END IF
  619.                   TJJ = CABS1( TJJS )
  620.                   IF( TJJ.GT.ONE ) THEN
  621. *
  622. *                       Divide by A(j,j) when scaling x if A(j,j) > 1.
  623. *
  624.                      REC = MIN( ONE, REC*TJJ )
  625.                      USCAL = ZLADIV( USCAL, TJJS )
  626.                   END IF
  627.                   IF( REC.LT.ONE ) THEN
  628.                      CALL ZDSCAL( N, REC, X, 1 )
  629.                      SCALE = SCALE*REC
  630.                      XMAX = XMAX*REC
  631.                   END IF
  632.                END IF
  633. *
  634.                CSUMJ = ZERO
  635.                IF( USCAL.EQ.DCMPLX( ONE ) ) THEN
  636. *
  637. *                 If the scaling needed for A in the dot product is 1,
  638. *                 call ZDOTU to perform the dot product.
  639. *
  640.                   IF( UPPER ) THEN
  641.                      CSUMJ = ZDOTU( J-1, A( 1, J ), 1, X, 1 )
  642.                   ELSE IF( J.LT.N ) THEN
  643.                      CSUMJ = ZDOTU( N-J, A( J+1, J ), 1, X( J+1 ), 1 )
  644.                   END IF
  645.                ELSE
  646. *
  647. *                 Otherwise, use in-line code for the dot product.
  648. *
  649.                   IF( UPPER ) THEN
  650.                      DO 130 I = 1, J - 1
  651.                         CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I )
  652.   130                CONTINUE
  653.                   ELSE IF( J.LT.N ) THEN
  654.                      DO 140 I = J + 1, N
  655.                         CSUMJ = CSUMJ + ( A( I, J )*USCAL )*X( I )
  656.   140                CONTINUE
  657.                   END IF
  658.                END IF
  659. *
  660.                IF( USCAL.EQ.DCMPLX( TSCAL ) ) THEN
  661. *
  662. *                 Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
  663. *                 was not used to scale the dotproduct.
  664. *
  665.                   X( J ) = X( J ) - CSUMJ
  666.                   XJ = CABS1( X( J ) )
  667.                   IF( NOUNIT ) THEN
  668.                      TJJS = A( J, J )*TSCAL
  669.                   ELSE
  670.                      TJJS = TSCAL
  671.                      IF( TSCAL.EQ.ONE )
  672.      $                  GO TO 160
  673.                   END IF
  674. *
  675. *                    Compute x(j) = x(j) / A(j,j), scaling if necessary.
  676. *
  677.                   TJJ = CABS1( TJJS )
  678.                   IF( TJJ.GT.SMLNUM ) THEN
  679. *
  680. *                       abs(A(j,j)) > SMLNUM:
  681. *
  682.                      IF( TJJ.LT.ONE ) THEN
  683.                         IF( XJ.GT.TJJ*BIGNUM ) THEN
  684. *
  685. *                             Scale X by 1/abs(x(j)).
  686. *
  687.                            REC = ONE / XJ
  688.                            CALL ZDSCAL( N, REC, X, 1 )
  689.                            SCALE = SCALE*REC
  690.                            XMAX = XMAX*REC
  691.                         END IF
  692.                      END IF
  693.                      X( J ) = ZLADIV( X( J ), TJJS )
  694.                   ELSE IF( TJJ.GT.ZERO ) THEN
  695. *
  696. *                       0 < abs(A(j,j)) <= SMLNUM:
  697. *
  698.                      IF( XJ.GT.TJJ*BIGNUM ) THEN
  699. *
  700. *                          Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
  701. *
  702.                         REC = ( TJJ*BIGNUM ) / XJ
  703.                         CALL ZDSCAL( N, REC, X, 1 )
  704.                         SCALE = SCALE*REC
  705.                         XMAX = XMAX*REC
  706.                      END IF
  707.                      X( J ) = ZLADIV( X( J ), TJJS )
  708.                   ELSE
  709. *
  710. *                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
  711. *                       scale = 0 and compute a solution to A**T *x = 0.
  712. *
  713.                      DO 150 I = 1, N
  714.                         X( I ) = ZERO
  715.   150                CONTINUE
  716.                      X( J ) = ONE
  717.                      SCALE = ZERO
  718.                      XMAX = ZERO
  719.                   END IF
  720.   160             CONTINUE
  721.                ELSE
  722. *
  723. *                 Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
  724. *                 product has already been divided by 1/A(j,j).
  725. *
  726.                   X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ
  727.                END IF
  728.                XMAX = MAX( XMAX, CABS1( X( J ) ) )
  729.   170       CONTINUE
  730. *
  731.          ELSE
  732. *
  733. *           Solve A**H * x = b
  734. *
  735.             DO 220 J = JFIRST, JLAST, JINC
  736. *
  737. *              Compute x(j) = b(j) - sum A(k,j)*x(k).
  738. *                                    k<>j
  739. *
  740.                XJ = CABS1( X( J ) )
  741.                USCAL = TSCAL
  742.                REC = ONE / MAX( XMAX, ONE )
  743.                IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN
  744. *
  745. *                 If x(j) could overflow, scale x by 1/(2*XMAX).
  746. *
  747.                   REC = REC*HALF
  748.                   IF( NOUNIT ) THEN
  749.                      TJJS = DCONJG( A( J, J ) )*TSCAL
  750.                   ELSE
  751.                      TJJS = TSCAL
  752.                   END IF
  753.                   TJJ = CABS1( TJJS )
  754.                   IF( TJJ.GT.ONE ) THEN
  755. *
  756. *                       Divide by A(j,j) when scaling x if A(j,j) > 1.
  757. *
  758.                      REC = MIN( ONE, REC*TJJ )
  759.                      USCAL = ZLADIV( USCAL, TJJS )
  760.                   END IF
  761.                   IF( REC.LT.ONE ) THEN
  762.                      CALL ZDSCAL( N, REC, X, 1 )
  763.                      SCALE = SCALE*REC
  764.                      XMAX = XMAX*REC
  765.                   END IF
  766.                END IF
  767. *
  768.                CSUMJ = ZERO
  769.                IF( USCAL.EQ.DCMPLX( ONE ) ) THEN
  770. *
  771. *                 If the scaling needed for A in the dot product is 1,
  772. *                 call ZDOTC to perform the dot product.
  773. *
  774.                   IF( UPPER ) THEN
  775.                      CSUMJ = ZDOTC( J-1, A( 1, J ), 1, X, 1 )
  776.                   ELSE IF( J.LT.N ) THEN
  777.                      CSUMJ = ZDOTC( N-J, A( J+1, J ), 1, X( J+1 ), 1 )
  778.                   END IF
  779.                ELSE
  780. *
  781. *                 Otherwise, use in-line code for the dot product.
  782. *
  783.                   IF( UPPER ) THEN
  784.                      DO 180 I = 1, J - 1
  785.                         CSUMJ = CSUMJ + ( DCONJG( A( I, J ) )*USCAL )*
  786.      $                          X( I )
  787.   180                CONTINUE
  788.                   ELSE IF( J.LT.N ) THEN
  789.                      DO 190 I = J + 1, N
  790.                         CSUMJ = CSUMJ + ( DCONJG( A( I, J ) )*USCAL )*
  791.      $                          X( I )
  792.   190                CONTINUE
  793.                   END IF
  794.                END IF
  795. *
  796.                IF( USCAL.EQ.DCMPLX( TSCAL ) ) THEN
  797. *
  798. *                 Compute x(j) := ( x(j) - CSUMJ ) / A(j,j) if 1/A(j,j)
  799. *                 was not used to scale the dotproduct.
  800. *
  801.                   X( J ) = X( J ) - CSUMJ
  802.                   XJ = CABS1( X( J ) )
  803.                   IF( NOUNIT ) THEN
  804.                      TJJS = DCONJG( A( J, J ) )*TSCAL
  805.                   ELSE
  806.                      TJJS = TSCAL
  807.                      IF( TSCAL.EQ.ONE )
  808.      $                  GO TO 210
  809.                   END IF
  810. *
  811. *                    Compute x(j) = x(j) / A(j,j), scaling if necessary.
  812. *
  813.                   TJJ = CABS1( TJJS )
  814.                   IF( TJJ.GT.SMLNUM ) THEN
  815. *
  816. *                       abs(A(j,j)) > SMLNUM:
  817. *
  818.                      IF( TJJ.LT.ONE ) THEN
  819.                         IF( XJ.GT.TJJ*BIGNUM ) THEN
  820. *
  821. *                             Scale X by 1/abs(x(j)).
  822. *
  823.                            REC = ONE / XJ
  824.                            CALL ZDSCAL( N, REC, X, 1 )
  825.                            SCALE = SCALE*REC
  826.                            XMAX = XMAX*REC
  827.                         END IF
  828.                      END IF
  829.                      X( J ) = ZLADIV( X( J ), TJJS )
  830.                   ELSE IF( TJJ.GT.ZERO ) THEN
  831. *
  832. *                       0 < abs(A(j,j)) <= SMLNUM:
  833. *
  834.                      IF( XJ.GT.TJJ*BIGNUM ) THEN
  835. *
  836. *                          Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
  837. *
  838.                         REC = ( TJJ*BIGNUM ) / XJ
  839.                         CALL ZDSCAL( N, REC, X, 1 )
  840.                         SCALE = SCALE*REC
  841.                         XMAX = XMAX*REC
  842.                      END IF
  843.                      X( J ) = ZLADIV( X( J ), TJJS )
  844.                   ELSE
  845. *
  846. *                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
  847. *                       scale = 0 and compute a solution to A**H *x = 0.
  848. *
  849.                      DO 200 I = 1, N
  850.                         X( I ) = ZERO
  851.   200                CONTINUE
  852.                      X( J ) = ONE
  853.                      SCALE = ZERO
  854.                      XMAX = ZERO
  855.                   END IF
  856.   210             CONTINUE
  857.                ELSE
  858. *
  859. *                 Compute x(j) := x(j) / A(j,j) - CSUMJ if the dot
  860. *                 product has already been divided by 1/A(j,j).
  861. *
  862.                   X( J ) = ZLADIV( X( J ), TJJS ) - CSUMJ
  863.                END IF
  864.                XMAX = MAX( XMAX, CABS1( X( J ) ) )
  865.   220       CONTINUE
  866.          END IF
  867.          SCALE = SCALE / TSCAL
  868.       END IF
  869. *
  870. *     Scale the column norms by 1/TSCAL for return.
  871. *
  872.       IF( TSCAL.NE.ONE ) THEN
  873.          CALL DSCAL( N, ONE / TSCAL, CNORM, 1 )
  874.       END IF
  875. *
  876.       RETURN
  877. *
  878. *     End of ZLATRS
  879. *
  880.       END
  881.